Digital equalizing method and radio communication receiver implementing said method

ABSTRACT

For processing samples of a signal received via a channel represented by an impulse response of W−1 coefficients, the method comprises: determining the W roots of the Z-transform of the channel impulse response; producing an intermediate signal by equalizing the received signal by a zero-forcing method or the like based on an impulse response whose Z-transform is a Z −1  polynomial of degree W-p having as roots those of the W roots which are furthest from the unit circle of the complex plane; and then obtaining estimations of the transmitted signal symbols by applying a Viterbi-type equalization method or the like based on an impulse response whose Z-transform is a Z −1  polynomial of degree p having as roots those of the W roots which are nearest to the unit circle.

BACKGROUND OF THE INVENTION

The invention concerns the digital equalization of signals. It has important applications in the field of radio communications.

The method is applied to the processing of a signal received from a transmitter over a channel between transmitter and receiver, whose response is known or has been estimated beforehand. One of the main issues is then that of the compromise between the performance of the equalizer and its complexity.

A complete maximum likelihood estimation of the transmitted signal is possible, for example by employing the Viterbi algorithm (see G. D. Forney Jr.: The Viterbi Algorithm, Proc. of the IEEE, Vol. 61, No. 3, March 1973, pages 268-278). Nevertheless, as soon as the impulse response of the channels becomes long or the number of possible discrete values becomes large, the exponentially increasing complexity of these methods renders them impracticable.

We consider the case of a radio communications channel used for the transmission of a signal composed of successive sequences or frames of n symbols, d_(k) (1≦k≦n). The symbols d_(k) have discrete values: binary (±1) in the case of a BPSK (binary phase shift keying) modulation; quaternary (±1+j) in the case of a QPSK (quaternary phase shift keying) modulation.

After baseband conversion, digitization and matched filtering, a vector Y of the received signal, corresponding to the symbols transmitted over the duration of a frame, is defined by the expression $\begin{matrix} {Y = {\left( \quad\begin{matrix} y_{1} \\ y_{2} \\ \quad \\ \vdots \\ \quad \\ y_{k} \\ \quad \\ \vdots \\ \quad \\ y_{L} \end{matrix}\quad \right) = {{{\left( \quad\begin{matrix} r_{0} & 0 & 0 & \cdots & 0 \\ r_{1} & r_{0} & 0 & \quad & \vdots \\ \quad & r_{1} & r_{0} & ⋰ & 0 \\ \vdots & \quad & r_{1} & ⋰ & 0 \\ \quad & \vdots & \quad & ⋰ & r_{0} \\ r_{w} & \quad & \vdots & \quad & r_{1} \\ 0 & r_{w} & \quad & \quad & \quad \\ 0 & 0 & r_{w} & \quad & \vdots \\ \vdots & \quad & ⋰ & ⋰ & \quad \\ 0 & \cdots & 0 & 0 & r_{W} \end{matrix}\quad \right)\left( \quad\begin{matrix} d_{1} \\ d_{2} \\ \quad \\ \vdots \\ \quad \\ d_{k} \\ \quad \\ \vdots \\ \quad \\ d_{n} \end{matrix}\quad \right)} + \left( \quad\begin{matrix} y_{N,1} \\ y_{N,2} \\ \quad \\ \vdots \\ \quad \\ y_{N,k} \\ \quad \\ \vdots \\ \quad \\ y_{N,L} \end{matrix}\quad \right)} = {{A.D} + Y_{N}}}}} & (1) \end{matrix}$ where W+1 is the length in numbers of bits of the estimated impulse response of the channel, r=(r₀, r₁, . . . , r_(w)) is the estimated impulse response of the channel, the r_(q) being complex numbers such that r_(q)=0 for q<0 or q>W, Y_(k) is the k-th complex sample received with 1≦k≦L=n+W, and Y_(N) is a vector of dimension L composed of samples of additive noise y_(N,k). The estimated impulse response r takes into account the propagation channel, the signal shaping by the transmitter and the receiver filtering.

The matrix A of dimension L×n has a Toeplitz-type structure, meaning that if α_(i,j) denotes the term in the i-th row and in the j-th column of the matrix A, then α_(i+1,j+1)=α_(i,j) for 1≦i≦L−1 and 1≦i≦n−1. The terms of the matrix A are given by: α_(1,j)=0 for 1<j≦n (A therefore has only zeros above its principal diagonal); α_(i,1)=0 for W+1>i≦L (band-matrix structure); and α_(i,1)=r_(i−1), for 1≦i≦W+1.

The matrix equation (1) expresses the fact that the signal received Y is an observation, with an additive noise, of the convolution product between the channel impulse response and the transmitted symbols. This convolution product can also be expressed by its Z-transform: Y(Z)=R(Z).D(Z)+Y_(N)(Z)  (2) where D(Z), Y(Z), R(Z) and Y_(N)(Z) are the Z-transforms of the transmitted symbols, of the received signal, of the impulse response and of the noise, respectively: $\begin{matrix} {{D(Z)} = {\sum\limits_{k = 1}^{n}\quad{d_{k} \cdot Z^{- k}}}} & (3) \\ {{Y(Z)} = {\sum\limits_{k = 1}^{L}\quad{y_{k} \cdot Z^{- k}}}} & (4) \\ {{R(Z)} = {\sum\limits_{q = 0}^{W}\quad{r_{q} \cdot Z^{- q}}}} & (5) \end{matrix}$

A conventional solution to solve a system such as (1) is the so-called “zero forcing” method, by which we determine the vector {circumflex over (D)}_(ZF) of n continuous components which minimizes the quadratic error ε=∥A{circumflex over (D)}−Y∥². Subsequently, a discretization of the components of vector {circumflex over (D)}_(ZF) relating to each channel is performed, often through a channel decoder. The least-squares solution {circumflex over (D)}_(ZF) is given by {circumflex over (D)}_(ZF)=(A^(H)A)⁻¹A^(H)Y, where A^(H) denotes the conjugate transpose of matrix A. We then are faced with the problem of inverting the Hermitian positive matrix ALA. The inversion can be effected by various classical algorithms, either directly (method of Gauss, Cholesky, etc.) or by an iterative technique (Gauss-Seidel algorithm, gradient algorithm, etc.).

The estimation error D−{circumflex over (D)}_(ZF) is equal to (A^(H)A)⁻¹A^(H)Y_(N), hence the solution includes noise with a variance: σ² =E(∥D−{circumflex over (D)} _(ZF)∥²)=N ₀×Trace[(A ^(H) A)⁻¹]  (6) where N₀ is the noise power spectral density. It can be seen that noise enhancement occurs when the matrix A_(H)A is badly conditioned, i.e. when it has one or more eigenvalues close to zero.

This noise enhancement is the main drawback of the conventional solution methods. In practice, the cases where the matrix A^(H)A is badly conditioned are frequent, especially with multiple propagation paths.

A relatively simple means to partly remedy this drawback is known, by accepting residual interference in the solution, i.e. by not adopting the optimum least-squares solution, but the solution: D_(MMSE)=(A^(H)A+{circumflex over (N)}₀)^(−A) ^(H)Y, where {circumflex over (N)}₀ denotes an estimation of the noise spectral density, that the receiver must then calculate. This is known as the MMSE (minimum mean square error) method and it allows the estimation variance to be reduced relative to the zero-forcing method, but introduces a bias.

The zero forcing methods and the like amount to performing an inverse filtering on the signal received by a filter which models the transfer function 1/R(Z) calculated using a certain approximation (quadratic in the case of zero forcing). Where one or more roots of the polynomial R(Z) (equation (5)) are situated on the unit circle, the theoretical inverse filter presents singularities rendering it impossible to estimate by a satisfactory approximation. In the case of the quadratic approximation, this corresponds to a divergence in the error variance σ² when the matrix A^(H)A has an eigenvalue equal to zero (relation (6)).

This problem does not arise in methods such as the Viterbi algorithm which intrinsically take the discrete nature of the symbols into account, but which require a much higher calculating power for large-sized systems.

U.S. Pat. No. 4,701,936 discloses a channel equalizer using an all-pass filter determined with reference to the Z-transform of the estimated channel impulse response.

An object of the present invention is to propose an equalization method achieving a good compromise between the reliability of the estimations and the complexity of the equalizer.

Another object is to produce an equalizer requiring a reasonable calculating power but capable of processing, with performance comparable to that of Viterbi equalizers, signals whose symbols have a relatively high number of states and/or signals carried on a channel with a relatively log impulse response.

SUMMARY OF THE INVENTION

Accordingly, the invention proposes a digital equalization method for estimating discrete information symbols from digital samples of a signal received over a transmission channel represented by a finite impulse response of W+1 coefficients, W being an integer greater than 1. This method comprises the steps of:

-   -   determining the W roots in the complex plane of the Z-transform         of the impulse response;     -   distributing the W roots into a first set of W-p roots and a         second set of p roots, p being an integer greater than 0 and         smaller than W, the roots of the second set being closer to the         unit circle than those of the first set according to a         determined distance criterion in the complex plane;     -   obtaining an intermediate signal by applying a first         equalization method to the received signal based on a finite         impulse response whose Z-transform, formed by a polynomial of         degree W-p in Z⁻¹, has roots which are the W-p roots of the         first set; and     -   obtaining estimations of the discrete information symbols by         applying a second equalization method to the intermediate signal         based on a finite impulse response whose Z-transform, formed by         a polynomial of degree p in Z⁻¹, has roots which are the p roots         of the second set.

The “first equalisation method” will usually be chosen so as to treat the unknown symbols as continuous variables. This leads to a process similar to an inverse filtering whose transfer function would have a form approaching the expression 1/RS(Z), where R^(S)(Z) denotes the polynomial in Z⁻¹ of order W-p having roots that are the W-p roots furthest from the unit circle in the complex plane. In particular, it could be of the “zero-forcing” type. This process generates only a small enhancement of the noise, since the roots of the associated transfer function in Z are relatively far from the unit circle.

For the p roots which are closest to the unit circle, we adopt measures to avoid or to limit the incidence of the noise enhancement problem. A MMSE or similar method can be chosen as the “second equalization method”. However, this second method will advantageous take into account the discrete nature of the unknown symbols. In particular, it may rely on a trellis algorithm, such as the Viterbi algorithm, whose implementation is common in channel equalizers where the system size is not too big.

The implementation of the second equalization method is generally more complex than the first. In each particular case, the choice of the number p allows the best compromise between the reliability of the estimations, which favors higher values of p, and the complexity of the equalizer, which favors lower values of p.

Another aspect of the present invention relates to a radio communications receiver comprising

-   -   conversion means to produce digital samples from a radio signal         received over a transmission channel represented by a finite         impulse response of W+1 coefficients, w being an integer greater         than 1;     -   means for measuring the channel impulse response;     -   means for calculating the W roots in the complex plane of the         Z-transform of the impulse response;     -   means for distributing the W roots into a first set of W-p roots         and a second set of p roots, p being an integer greater than 0         and smaller than W, the roots of the second set being closer to         the unit circle than those of the first set according to a         determined distance criterion in the complex plane;     -   a first equalization stage for producing an intermediate signal         by applying a first equalization method to the received signal         based on a finite impulse response whose Z transform, formed by         a polynomial of degree W-p in Z⁻¹, has roots which are the W-p         roots of the first set; and     -   a second equalization stage for producing estimations of the         discrete symbols of a signal carried on the channel by applying         a second equalization method to the intermediate signal based on         a finite impulse response whose Z transform (R^(I)(Z)), formed         by a polynomial of degree p in Z⁻¹, has roots which are the p         roots of the second set.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram of an example of a radio communications receiver according to the invention;

FIG. 2 is a flow chart showing an embodiment of the method according to the invention; and

FIG. 3 is a diagram illustrating the performances

DESCRIPTION OF PREFERRED EMBODIMENTS

The receiver represented in FIG. 1 comprises a radio stage 1 which receives the signal picked up by antenna 2 and converts it to the baseband. The baseband signal is digitized by the analog-to-digital converter (ADC) 3 and fed to the reception filter 4. Filter 4 provides matched filtering with respect to the signal shaping by the transmitter. It delivers a digital signal at a rate of one complex sample per transmitted symbol.

The digital signal is fed to a demodulator comprising, on the one hand, a synchronization and channel estimation module 6 and, on the other hand, an equalizer 7.

The synchronization and channel estimation are, for example, performed conventionally by means of a synchronization sequence included by the transmitter in each signal frame. Detection of this sequence, known to the receiver, allows, on the one hand, the receiver to be synchronized relative to the temporal structure of the transmitted frames, and on the other to estimate the impulse response r=(r₀, r₁, . . . , r_(w)) of the channel carrying the transmitted frames. The impulse response calculated by module 6 is fed to the equalizer 7.

Equalizer 7 functions for example according to the flow chart shown in FIG. 2 to process each frame of the received signal, which are presented in the form of a vector ${Y = \left( \quad\begin{matrix} y_{1} \\ \vdots \\ y_{L} \end{matrix}\quad \right)},$ with L=n+W, adopting the previous notation.

The channel estimation module 6 having supplied the W+1 complex coefficients r_(q) of the estimated channel impulse response, the first step 10 consists of a search of the W roots of the Z-transform of the impulse response, given by equation (5). Various conventional methods for the determination of complex roots of a polynomial can be employed at step 10. To this end, reference can be made to the work of E. DURAND: “Numerical Solutions of Algebraic Equations, Volume I, Equations of the Type F(x)=0”, Editions Masson, 1960.

The W roots thus found α₁, α₂, . . . , α_(W) are then ordered to allow distribution into two sets, one containing the W-p roots furthest from the unit circle, and the other the p roots closest to the unit circle.

To this end, a distance δ_(q) is calculated at step 11 for each of the roots α_(q) (1≦q≦W). This distance can be conveniently obtained as follows $\begin{matrix} {\delta_{q} = \left\{ \quad\begin{matrix} {1 - {\alpha_{q}}} & {{{si}{\alpha_{q}}} \leq 1} \\ {1 - {1/{\alpha_{q}}}} & {{{si}{\alpha_{q}}} > 1} \end{matrix}\quad \right.} & (7) \end{matrix}$

In step 12, the roots aq of the transfer function R(Z) are sorted in order of decreasing distance: δ₁≧δ₂ . . . ≧δ_(W). Then, the first W-p roots which are the furthest from the unit circle, are separated from the remaining roots α_(w-p+1), . . . , α_(W).

In step 13, the equalizer 7 develops a polynomial in Z⁻¹ defined by: $\begin{matrix} {{R^{S}(Z)} = {{\prod\limits_{q = 1}^{W - p}\quad\left( {1 - {\alpha_{q} \cdot Z^{- 1}}} \right)} = {\sum\limits_{q = 0}^{W - p}\quad{s_{q} \cdot Z^{- q}}}}} & (8) \end{matrix}$

This allows the coefficients s_(q) of the transfer function R^(S)(Z) associated with the impulse response s=(s₀, s₁, . . . , s_(W-p)) of a virtual channel to be determined, which corresponds to the estimated transmission channel with the contributions closest to the singularity zones removed.

We then proceed to a first equalization 14 which is effectively an inverse filtering approaching the transfer function 1/R^(S)(Z). Several implementations can be used to perform this inverse filtering. In particular, a zero-forcing equalization can be performed, as previously discussed. Regarding these methods, reference can be made to the work of J. G. Proakis: “Digital Communications”, McGraw-Hill, 2^(nd) Edition, 1989.

The inverse filtering 14 produces an intermediate signal in the form of a vector Y′ of L′=n+p samples y′₁, . . . , y′_(L). In the case of the zero-forcing method, the vector Y′ is obtained from the matrix equation: Y′(A′ ^(H) A′)⁻¹ A′^(H) Y  (9)

In expression (9), A′ denotes a matrix having a Toeplitz structure with n+W rows and n+p columns, generated from the coefficients s_(q) of the polynomial R^(S)(Z): $\begin{matrix} {A^{\prime} = \left( \quad\begin{matrix} s_{0} & 0 & 0 & \cdots & 0 \\ s_{1} & s_{0} & 0 & \quad & \vdots \\ \quad & s_{1} & s_{0} & ⋰ & 0 \\ \vdots & \quad & s_{1} & ⋰ & 0 \\ \quad & \vdots & \quad & ⋰ & s_{0} \\ s_{w - p} & \quad & \vdots & \quad & s_{1} \\ 0 & s_{w - p} & \quad & \quad & \quad \\ 0 & 0 & s_{w - p} & \quad & \vdots \\ \vdots & \quad & ⋰ & ⋰ & \quad \\ 0 & \cdots & 0 & 0 & s_{w - p} \end{matrix}\quad \right)} & (10) \end{matrix}$

The sorting of the roots am results in the eigenvalues of matrix A′^(H)A′ being relatively far removed from 0.

Alternatively, the inverse filtering can be achieved by cascading W-p filtering cells, each corresponding to the inverse of a transfer function R_(q) ^(s)(Z)=1−α_(q)Z⁻¹, for 1≦q≦W-p. If |αq|=1, the inverse filter of R_(q) ^(s)(Z) is unfeasible. If |α_(q)|<1, we can develop 1/R_(q) ^(s)(Z) in the form: $\begin{matrix} {\frac{1}{R_{q}^{S}(Z)} = {1 + {\alpha_{q} \cdot Z^{- 1}} + {\alpha_{q}^{2} \cdot Z^{- 2}} + \ldots + {\alpha_{q}^{m} \cdot Z^{- m}} + \ldots}} & (11) \end{matrix}$ Development (11) is causal and stable since the convergence domain contains the unit circle. The inverse filter cell can thus be generated in either a transverse or a recursive form.

If ⊕α_(q)|>1, then 1/R_(q) ^(s)(Z) can be developed in the form: $\begin{matrix} {\frac{1}{R_{q}^{S}(Z)} = {{- \alpha_{q}^{- 1}} \cdot Z \cdot \left( {1 + {\alpha_{q}^{- 1} \cdot Z} + {\alpha_{q}^{- 2} \cdot Z^{2}} + \ldots + {\alpha_{q}^{- m} \cdot Z^{m}} + \ldots} \right)}} & (12) \end{matrix}$ Development (12) is anti-causal and stable. To perform the inverse filtering cell, the development (12) is truncated and an transverse filter implementation is adopted. The anti-causality results in a delay corresponding to the length of the retained response.

We note that developments (11) and (12) justify the unit circle distance criterion δ_(q) used according to relation (7).

In step 15, the equalizer 7 develops a polynomial of degree p in Z⁻¹, whose roots correspond to the p roots of R(Z) closest to the unit circle, such that R(Z)=R^(S)(Z).R^(I)(Z): $\begin{matrix} {{R^{I}(Z)} = {{r_{0} \cdot {\prod\limits_{q = {W - p + 1}}^{W}\quad\left( {1 - {\alpha_{q} \cdot Z^{- 1}}} \right)}} = {\sum\limits_{q = 0}^{p}\quad{t_{q} \cdot Z^{- q}}}}} & (13) \end{matrix}$

The complex coefficients t_(q) define the impulse response of another virtual transmission channel, whose equalization by a zero-forcing method or the like would pose problems of noise enhancement.

The intermediate signal Y′ is then subjected to an equalization according to another method, based on the impulse response t=(t⁰, t₁ . . . , t_(p)). This second equalization 16 can be conveniently performed by means of a Viterbi trellis (see above-cited article from G. D. Forney Jr., or the above-cited work of J. G. Proakis).

The second equalization stage 16 generates estimations d_(k) of the symbols of the frame (1≦k≦n). These estimations d_(k) constituting the output of equalizer 7 can be input to a deinterleaving module 8 then to a channel decoder 9 which detects and/or corrects possible transmission errors.

FIG. 3 illustrates the performance of the method in the case of the transmission of a signal frame according to the GSM European mobile telephone system, replacing the binary modulation of the GMSK type by an eight state phase modulation (8-PSK modulation). The impulse response was truncated at five bit time periods (W=4). FIG. 3 shows the dependence of the bit error rate (BER), expressed in %, on the ratio E_(b)/N₀ between the energy-per-bit and the spectral density of the noise, expressed in decibels (dB). The BER is that observed in the symbol estimations after deinterleaving and channel decoding carried out by a method conforming to those employed in the GSM standard. Curve I shows the results obtained by the pure zero-forcing method, i.e. in the limit case where p=0. Curve II shows the theoretical result that would be obtained by equalizing the channel purely with the Viterbi algorithm (limit case where p=W). In practice, the corresponding trellis should have 8⁴=4096 states, meaning that the Viterbi equalizer could not be realized using current technology. The difference between curves I and II illustrates the superiority of the Viterbi algorithm which produces maximum likelihood estimations.

Curves III and IV show the results obtained by the method according to the invention, in the cases where p=1 and p=2, respectively. The quite significant improvement already achieved for the value p=1 with respect to pure zero forcing can be seen.

By way of an example, equalization of a GSM signal frame by a pure Viterbi algorithm, under the conditions of FIG. 3, would require of the order of 8.45 million floating-point operations, i.e. around 1.83 Gflops, while the implementation of the present invention under the same conditions requires of the order of 19,000 floating-point operations (≈4.2 Mflops) in the case of p=1, including searching the roots of R(Z) and 1/RS(Z) inverse filtering by the zero-forcing method. The number increases to the order of 129,000 operations (≈28 Mflops) in the case of p=2, which remains compatible with the capacity of currently available digital signal processors (DSP). 

1. A digital equalization method for estimating discrete information symbols from digital samples of a signal received over a transmission channel represented by a finite impulse response of W+1 coefficients, W being an integer greater than 1, comprising the steps of: determining W roots in the complex plane of a Z-transform of the impulse response; distributing the W roots into a first set of W-p roots and a second set of p roots, p being an integer greater than 0 and smaller than W, the roots of the second set being closer to a unit circle of the complex plane than those of the first set according to a determined distance criterion in the complex plane; obtaining an intermediate signal by applying a first equalization method to the received signal based on a finite impulse response having a Z-transform consisting of a polynomial of degree W-p in Z⁻¹, having roots equal to the W-p roots of the first set; and obtaining estimations of the discrete information symbols by applying a second equalization method to the intermediate signal based on a finite impulse response having a Z-transform consisting of a polynomial of degree p in Z⁻¹, having roots equal to the p roots of the second set.
 2. A method according to claim 1, wherein the first equalization method yields the intermediate signal in the form of a vector Y′ of n+p samples obtained according to the relation: Y′=(A′ ^(H) A′)⁻¹ A′ ^(H) Y where n is an integer representing a number of the discrete information symbols, Y is a vector composed of n+W samples of the received signal, and A′ is a matrix with n+W rows and n+p columns having a Toeplitz structure formed from the coefficients of said polynomial of degree W-p in Z⁻¹, and H denotes conjugate transpose.
 3. A method according to claim 1, wherein the second equalization method comprises implementing a Viterbi algorithm.
 4. A method according to claim 1, wherein the unit circle distance criterion, used to distribute the W roots α₁, . . . , α_(W) of the Z-transform of the channel impulse response into the first and second sets, is expressed as a distance δ_(q) of the form δ_(q)1−|α_(q)| if |α_(q)|≦1, and of the form δ_(q)=1−1/|α_(q)| if |α_(q)|>1, for 1≦q≦W.
 5. A radio communications receiver comprising: conversion means to produce digital samples from a radio signal received over a transmission channel represented by a finite impulse response of W+1 coefficients, W being an integer greater than 1; means for measuring the channel impulse response; means for calculating W roots in a complex plane of the Z-transform of the impulse response; means for distributing the W roots into a first set of W-p roots and a second set of p roots, p being an integer greater than 0 and smaller than W, the roots of the second set being closer to a unit circle of the complex plane than those of the first set according to a determined distance criterion in the complex plane; a first equalization stage for producing an intermediate signal by applying a first equalization method to the received signal based on a finite impulse response having a Z-transform consisting of a polynomial of degree W-p in Z⁻¹, having roots equal to the W-p roots of the first set; and a second equalization stage for producing estimations of discrete symbols of a signal carried on the channel by applying a second equalization method to the intermediate signal based on a finite impulse response having a Z-transform consisting of a polynomial of degree p in Z⁻¹, having roots equal to the p roots of the second set.
 6. A receiver according to claim 5, wherein the first equalization stage is arranged to yield the intermediate signal in the form of a vector Y′ of n+p samples obtained according to the relation: Y′=(A′ ^(H) A′)⁻¹ A′ ^(H) Y where n is an integer representing a number of the discrete symbols, Y is a vector composed of n+W samples of the received signal, and A′ is a matrix with n+W rows and n+p columns having a Toeplitz structure formed from the coefficients of said polynomial of degree W-p in Z⁻¹, and H denotes conjugate transpose.
 7. A receiver according to claim 5, wherein the second equalization stage is arranged to implement a Viterbi algorithm.
 8. A receiver according to claim 5, wherein the means for distributing the W roots α₁, . . . , α_(W) into the first and second sets make use of a unit circle distance criterion expressed as a distance δ_(q) of the form δ_(q)=1−|α_(q)| if |α_(q)|≦1, and of the form δ_(q)=1−1/|α_(q)|if |α_(q)|>1, for 1≦q≦W. 